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Abstract 


Full seven-equation Reynolds stress turbulence models are promising tools for today’s aerospace 
technology challenges. This paper examines two such models for computing challenging turbulent flows 
including shock-wave boundary layer interactions, separation and mixing layers. The Wilcox and the 
SSG/LRR full second-moment Reynolds stress models have been implemented into the FUN3D 
unstructured Navier-Stokes code and were evaluated for four problems: a transonic two-dimensional 
diffuser, a supersonic axisymmetric compression corner, a compressible planar shear layer, and a 
subsonic axisymmetric jet. Simulation results are compared with experimental data and results computed 
using the more commonly used Spalart-Allmaras (SA) one-equation and the Menter Shear Stress 
Transport (SST-V) two-equation turbulence models. 


Nomenclature 


aij Reynolds-stress anisotropy tensor 

b mixing layer thickness 

Cr skin friction coefficient 

»C; coefficients of the Speziale-Sarkar-Gatski model 
i coefficients of the Wilcox RSM pressure-strain 
Cy, equilibrium parameter 

D Reynolds stress diffusion coefficient 

Dij specific diffusion tensor 

d distance to nearest wall point 

Fy Menter’s blending function 

jel local diffuser height 

Atiroa diffuser throat height 

k specific turbulent kinetic energy 

Mi; turbulent mass flux tensor 

P, specific turbulence production tensor 

static pressure 

radial coordinate 


P 

r 

Sij strain rate tensor 

Sij traceless strain rate tensor 
t 


time 
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velocity component in the x-direction 
local mean streamwise velocity 


u 

U 

At freestream velocity difference, A“ = i, — tz 

Uj velocity component 

Wij rotation tensor 

x spatial coordinate in the streamwise or axial direction 

Xi Cartesian coordinates 

y spatial coordinate in the vertical or wall-normal direction 

Vo mixing layer centerline (midpoint between locations where t% = ti, — 0.1At% and 
i = ti, + 0.1Ai) 


a coefficient of @ production 

a closure coefficient 

Ay coefficient of @ production 

B coefficient of m destruction 

B closure coefficient 

Bo closure coefficient 

Bo coefficient of @ destruction 

V closure coefficient 

Oij Kronecker delta 

Eij turbulent dissipation rate 

¢ argument to F, blending function 

m dynamic viscosity 

Ij; specific pressure-strain correlation tensor 

p density 

o coefficient of diffusion 

o* closure coefficient 

Oa coefficient of cross-diffusion 

Ow coefficient of diffusion 

Tij viscous stress tensor 

© coefficient of € equation 

go coefficient of w equation 

w specific dissipation rate 

Accents 

gp" fluctuating component, 6” = ¢-— 

z z . 1 ptotAt 

p Reynolds averaged component, @ = lim — pdt 
At At “to 

x zy 14... 1 ttt 

p Favre averaged component, @ = 5 jim, re me podt 

Subscripts 

T turbulent 

t total 

0 stagnation condition 

oo freestream condition 
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1.0 Introduction 


For aerospace propulsion flows, the most common practice in computational fluid dynamics (CFD) 
analyses is to use Reynolds-Averaged Navier-Stokes (RANS) codes with one- or two-equation turbulence 
models. Current RANS turbulence models predict steady, fully turbulent attached flows at all speed 
regimes reasonably well, but are still unable to reliably predict most separated flows. Large-eddy 
simulation (LES) and direct-numerical simulation (DNS) methods are being used for some applications, 
however they require very fine grids for wall bounded flows and at high Reynolds numbers, and therefore 
will not be practical for many years (Refs. 1 and 2). Hybrid unsteady RANS/LES methods are 
increasingly common for certain classes of simulations including separated flows, although the techniques 
to combine the near-wall RANS region with the outer, large-eddy simulation region need further 
development (Ref. 1). 

For aerospace propulsions flows, RANS will be used for a significant portion of CFD analyses for the 
foreseeable future, due to limitations in computational power (Ref. 1). Traditional RANS linear and 
nonlinear one- and two-equation turbulence models, such as Spalart-Allmaras (SA) (Ref. 3) and Menter 
Shear-Stress Transport (SST) (Ref. 4), are frequently used. These linear models incorporate the 
Boussinesq approximation to give the Reynolds shear stress tensor in terms of the mean strain rate tensor 
and the eddy viscosity (Ref. 5). Nonlinear models incorporate additional higher-order terms that are 
functions of the mean strain and rotation rate tensors. An alternative to these approaches is to use a more 
advanced form of RANS turbulence modeling, where the individual transport equations for the Reynolds 
stresses are solved, allowing for a more detailed representation of the flow physics. These models are 
known as full second-moment Reynolds Stress Models (RSMs). 

This paper describes a study using the FUN3D (Ref. 6) unstructured CFD code to evaluate two 
RSMs: the Wilcox RSM (Ref. 5) and the combined Speziale-Sarkar-Gatski/Launder-Reece-Rodi 
(SSG/LRR) RSM (Ref. 7). This work supports the NASA Revolutionary Computational Aerosciences 
technical challenge to, “/dentify and down-select critical turbulence, transition, and numerical method 
technologies for 40% reduction in predictive error against standard test cases for turbulent separated 
flows, evolution of free shear flows and shock-boundary layer interactions on state-of-the-art high 
performance computing hardware.” (Ref. 1). The first step in this study was to reproduce previous results 
obtained with the Wind-US code using the SA and SST turbulence models in order to confirm that the 
turbulence models are implemented equivalently in each code. Next, FUN3D was run with the two 
RSMs. Four cases were examined: a transonic diffuser (Refs. 8 to 10), supersonic flow over an 
axisymmetric compression corner (Refs. 11 to 13), a compressible planar shear layer (Ref. 14), and a 
subsonic axisymmetric jet (Refs. 11 and 15). Results were compared with the SA and SST models, and 
with experimental data. 


2.0 The Codes 


The FUN3D code was used for the computations described herein and compared with 
results previously computed using the Wind-US code. Both are NASA production codes using 
Reynolds-averaged flow solvers and are described briefly below. 


2.1 Wind-US 


The Wind-US code (Refs. 16 and 17) has both structured and unstructured solvers. The flow 
equations are evaluated using second-order-accurate finite differences defaulting to second-order Roe for 
structured grids and second order HLLE for unstructured grids, although other schemes may be specified 
in the user inputs. The partial differential equations are written in conservative form and explicit terms are 
modeled using either upwind or central differencing. Implicit terms are computed using either an 
approximately factored or four-stage Runge-Kutta scheme. 
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The Wind-US results described in this paper are previously reported computations using the 
well-known Spalart-Allmaras (SA) (Ref. 3) one-equation and Menter shear-stress transport with vorticity 
source term (SST-V) (Ref. 4) two-equation turbulence models. The results are compared with FUN3D 
results computed using the same turbulence models in order to verify that these baseline turbulence 
models are implemented equivalently in each code. Additional details can be found at 
www.erc.nasa.gov/www/winddocs and in Reference 18. 


2.2. FUN3D 


FUN3D (Ref. 6) is a node-centered, unstructured implicit solver developed by researchers at the 
NASA Langley Research Center. It uses finite-volume discretization and is formally second-order 
accurate in space. Explicit terms are calculated using Roe’s flux difference splitting, however other 
methods are available. More information about FUN3D can be found at fun3d.larc.nasa.gov and in 
Reference 6. Note that at the time of this study, the full Reynolds stress models in FUN3D cannot be used 
for periodic grids whose side planes are rotated through a small angle. Instead, it requires 90-degree grids 
with both axes aligned with one of the constant x-, y- and z-coordinate surfaces. For the SA and SST-V 
models, small-angle periodic boundary conditions are available, so grids one-cell wide in the 
circumferential direction were used. 


2.3 Solver Inputs 


The inputs to both codes were set at the values suggested for each case according to the code 
documentation or recommended practices. For Wind-US, the right-hand-side inviscid explicit operator 
was set to the second-order HLLE flux splitting algorithm for the transonic diffuser case which was 
computed using an unstructured grid and the second-order Roe upwind scheme for the other cases, which 
were computed using structured grids. For FUN3D, the Roe second order upwind scheme was used. All 
cases were run steady state, with the exception of the axisymmetric jet, which required a time-accurate 
solution to reach convergence. 


3.0 The Turbulence Models 


The full RSM turbulence models evaluated in this work are the combined Speziale-Sarkar-Gatski/ 
Launder-Reece-Rodi (SSG/LRR) RSM and the Wilcox RSM. The more standard Spalart-Allmaras (SA) 
one-equation model (Ref. 3) and the Menter shear stress transport (SST-V) two-equation model with a 
vorticity source term (Ref. 4) were also used for comparison with the RSMs. Details about all of these 
models are available on the Turbulence Modeling Resource (TMR) website, turbmodels.larc.nasa.gov 
(Ref. 19). For one case, the mixing layer flow, a nonlinear explicit algebraic stress model (EASM), available 
in Wind-US, was also used (Ref. 20). The SA and SST-V models use the Boussinesq approximation, 


7 ~ 1s 2 
Pi; = 2Ur (Sj — =5iu5ij) — 3 Pk; (1) 
a constitutive relation, to compute the turbulent stress tensor, T; rp defined below as, 


Tij at Sup ae (2) 


On the other hand, the RSMs solve transport equations for each of the six unique Reynolds stresses. A 
seventh equation is required to determine the length scale variable. 
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The SSG/LRR full Reynolds stress model, as described on the TMR, and in References 7 and 21, was 
developed under the European Union project FLOMANIA and consists of the six equations of the 
Reynolds stress transport equations, plus Menter’s baseline w equation for the length scale. The 
Reynolds-stress transport equation is given by Equation (3). 


O(ptiy) , A(ptyju 
ig ER —pPij — pli + peiy — PDijy — PMi; (3) 


The pressure-strain model, given by Equation (4), blends the Speziale-Sarkar-Gatski model with the 
Launder-Reece-Rodi model near the wall, as defined in References 7 and 19. 


1 * ~ ~ ~ dines ~ * ~ ~ Tox 
Nij = (ie + = Ci Pax) aij + C2E (Gide; = + G1 15ij) + (C3 = C3 Ger Gye KS; + 
~f{. = ee Quis ox ~p ee 
Cyk (Gadi + Gj Siz — = GS yj) + Csk(GinWjx + Gj Wir) (4) 


The remaining terms in Equation (3) are given below. 


Pij = Tix 3 Mar, Production (5) 
Peij = =Pdijé Dissipation (6) 
aij = — a - = 6; j Aniosotropy Tensor (7) 
S; j= . (= + =| Strain Rate Tensor (8) 
Sij = Si - = SKK 6ij Traceless Strain Rate Tensor (9) 
W; j= ~ (= = =) Averaged Rotation Tensor (10) 
BDi = — 5 (iu — p Putt) wal Diffusion (1) 


The turbulent mass flux term, M;;, is assumed to be negligible. 

A length scale equation is required to close the system. The SSG/LRR model uses an w-equation 
similar to Menter’s SST-V model which blends an w-equation near the wall with an ¢ equation in the 
outer boundary layer. 


A(Pw) , A(Pwitk) _ , WPPKk pp ~ 24 9 [(- ze) 0) B nay (2k aw ) 
ot + OXR 7 %w k 2 Bu po + Ox [(Z = Pw OXk + Oa (co) ae (= Ox,’ 0 (12) 
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The isotropic dissipation rate is defined as 


e=C,ko (13) 


with C,, = 0.09. The equations used for blending the coefficients @ = a4, Bw,w,Oq are given below. 


$= Fo + (1-Fy)o© (14) 
F, = tanh(¢*) (15) 
oe Vk  500u 40 pk 
gomnin as (saa Gott) Tp a oe 
dw OxXROX,’ 


where d is the distance to the nearest wall. The pressure-strain coefficients are blended between 
Launder-Reece-Rodi (LRR) (Ref. 22) near walls (without wall-correction terms) and Speziale-Sarkar- 
Gatski (SSG) (Ref. 23) away from walls. The coefficients are given in Table 1. 

The diffusion term given in Equation (11), above, is a generalized gradient diffusion model (Ref. 24). 
In the FUN3D input file, the SSG/LRR RSM is specified as “SSG/LRR-RSM-w2012.” There is also a 
version of the SSG/LRR RSM which uses a simple diffusion model (Ref. 25) and is specified as 
“SSG/LRR-RSM-w2012-SD” in the input file. In this model, the diffusion term is modeled as given in 
Equation (17) 


es (# -2 tr) al Simple Diffusion (17) 


with: 
D =0.5C,F, + =0.22(1 — Fy) (18) 


Both diffusion terms are used in the paper. The transonic diffuser and mixing layer simulations both used 
the standard diffusion model. The simple diffusion model, which tends to be more stable, was used for the 
shock-wave boundary-layer interaction and the axisymmetric jet, because these simulations had more 
difficulty reaching convergence. 

The Wilcox RSM as described on the TMR and in Reference 5 also solves the six equations for the 
Reynolds stress tensor and an omega equation for the specific dissipation rate, w. The Reynolds stress 
equation is given by Equation (19). 


ED 4 See) —pPyj — ally +2B*pwk6;,; + = [e+ 09 =] (19) 
TABLE 1.—COEFFICIENTS FOR THE SSG/LRR FULL REYNOLDS STRESS TURBULENCE MODEL 
ty Bu Ow Og Cy Cy C2 C3 C3 Cy Cs D 
LRR © | 0.5556 | 0.075 | 0.5 0 18 | 0 0 | 08 | 0 | (ckRR +6) (—7C#R10)| 0.75C, 
11 11 
SSG & 0.44 0.0828 0.856 1.712 1.7 0.9 1.05 0.8 0.65 0.625 0.2 0.22 
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The specific dissipation rate is given by Equation (20). 


Opw O(Pwit;) _ pw Ou; BD 2 Pp 0k Ow 


= a—T,—-— ee 
@ Ox; OX; 


) = Ow 
at ax; k 4 ax; % alG + One) al (20) 


The pressure strain correlation is given by: 


Tj = B* Co (ti +265) = a(Py = = P6i;) —B (Di — = P6i;) = pk (Si =a = Six ij) (21) 
with, 
i. _ Oil, di, 
PS See Lr = pk/w Dij = Tik Ox, + Tk 3, (22) 


The production, P;;, given by Equation (5), the simple strain rate tensor, 5;;, given by Equation (8), 
and the averaged rotation tensor, W; j» given by Equation (10), are used in both the SSG/LRR and Wilcox 
RSMs. The closure coefficients for Equation (21) are given in Table 2, and in Equations (23) and (24). 


© aj 9%) = 23 
Oise dk 80 . 4 (23) 
” Ox; Ox; 
1+85X,, WijW jxSki A 1 dtm 
= = [Wu seo et Sa S508 ay, 24 
fp 1+100X,5 Xu (B*w)3 ki ki 2 0xm k,i ( ) 


The two-equation EASM (Ref. 20), available in Wind-US, was used for comparison with the other 
models for the mixing layer case. This model is derived from a reduced form of the Reynolds stress 
transport equations. The resultant expression for the Reynolds stresses is similar to Equation (1), but 
includes terms that are non-linear in the strain and rotation rate tensors. This particular EASM uses a 
simplified form of the SSG pressure-strain relation that is a linear function of the anisotropy. EASMs 
provide more information about the shear stress than the SA and SST-V models, while remaining less 
computationally expensive than the RSMs. Wind-US EASM results are shown for the mixing layer case 
in Section 4.3. 


TABLE 2.—CLOSURE COEFFICIENTS FOR THE WILCOX FULL REYNOLDS-STRESS TURBULENCE MODEL 
a B 4 C, C2 a B B* o o* Bo 
(8+¢,)/11 | (8-c)/11_ | (0c, —4)/55 SRE et ope Bell Sel Fog | oee || aovOR 
5 19 25 100 
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4.0 Test Cases and Results 


4.1 Transonic Diffuser 


4.1.1 Experimental Setup 


The first case examined was transonic flow through a converging-diverging diffuser. This flow has a 
strong normal shock wave in the throat causing a separated region to form on the top wall. This case was 
run previously with Wind-US, and was used to verify that the SA and SST-V turbulence models produce 
the same result in both Wind-US and FUN3D. Once this verification was complete, the flow was 
computed using the RSMs. Comparisons were made with the experimental data obtained by Sajben et al. 
(Ref. 8). The duct geometry in the x-y plane is shown along with Mach contours computed using the SA 
model in Figure 1. The geometry is a convergent/divergent channel with a flat bottom and a contoured top 
wall. The flow conditions at the inflow plane are Mach 0.9 with a total pressure and temperature of 
19.58 psi and 540 °R. A normal shock wave forms downstream of the throat at x/Hrprogt = 1.98, 
causing the flow on the top wall to separate. It reattaches at approximately six throat heights downstream. 


4.1.2 CFD Methodology 


The unstructured grid used for this case, shown in Figure 2, contains hexahedral cells near the walls 
and tetrahedral cells in the center. A constant area section extends 10 throat-heights downstream of the 
diffuser exit station to eliminate any boundary layer separation at the outflow as the flow sets up. (Only a 
portion of this is shown in Figure 2.) The grid contains 54,854 points, with 301 points in the axial 
direction. In the vertical direction, there are 41 vertical grid points in the hexahedral grid near each wall 
and approximately 91 total points in the vertical direction. The grid is one throat height deep in the 
transverse direction. The y" values one point from the wall are approximately 1.0. This grid is also 
described on the NPARC Alliance Verification and Validation Website (Ref. 26). Results computed on 
this grid were nearly equivalent to those computed on a grid twice as dense. 

The boundary conditions used are those described in Reference 26 to set up a strong shock in the 
diffuser. The top and bottom walls are adiabatic and viscous. The inflow conditions are set to the 
Mach 0.9 conditions described above. The outflow static pressure is set to 14.775 psi to allow the normal 
shock to form at the desired location in the throat, x/Htnrog¢ = 1.98. The full diffusion version of the 
SSG/LRR RSM was used. 


2 Mach 
number 
- ° 1.3 
8 12 
24 1.0 
z 0.8 
> 0.6 
0 0.4 
0.2 
4 —2 0 2 4 6 8 
X/Hthroat 


Figure 1.—Transonic diffuser. Mach contours computed with Wind-US and the SA turbulence model. 


y/Hthroat 


X/Hthroat 


Figure 2.—Transonic diffuser. Computational mesh. 
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4.1.3 Results 


The results computed using Wind-US and FUN3D with the SA and SST-V models are shown in 
Figure 3 and Figure 4. The normalized pressure on the top and bottom walls in Figure 3 indicates that 
both codes give nearly equivalent results. Both codes and turbulence models predict the pressure rise at 
the normal shock wave. However, just downstream of the shock, at x/Htnyoqge between 3 and 6, all of the 
results slightly overshoot the experimental values. Further downstream, they agree with the experiment. 
The velocity profiles shown in Figure 4 also indicate nearly equivalent results between the two codes. At 
X/Henroat = 2-9 and 4.6, the CFD results are all correctly predicting that the flow is separated on the top 
wall, and the velocity profiles are matching up well with the experiment. At x/Htnroat = 6.4, however, 
the CFD results still indicate that the flow on the top wall is separated whereas the experiment indicates 
that it is attached. At x/Htnroat = 7.5, the flow is attached for both the CFD and experiment, however, 


the CFD results predict thinner boundary layers on the top and bottom walls. 


3 =-1 


X/Athroat throat = 0.1445 ft 


1 3 5 


o 1.0 ° Experiment g 18 
Za 0.9 | Bottom wall Wind-US SA a= 09 Top wall 
RO — - + Wind-US SST-V ee” 
fo 08! | FUN3D SA fio 08 ¥ 
25 — — FUN3D SST-V 2s 
ES a 0.7 g 7 0.7 
(e) 
ice ¥ Sos 
sa ¥: o2Q ~ 
ES 04] Eo 04 
ro) Oo 
Zz 03 z 0.3 


5 3 1 1 3 5 7 @Q 
X/Athroat throat = 0.1445 ft 


Figure 3.—Transonic diffuser. Normalized wall pressure on the bottom and top walls computed with Wind-US and 


FUNSD. 

1.2 X/H = 2.9 1.2/5 7 x/H = 4.6 
= 1.0 — 1.0 
SE og SE og 
2D * 2D 
iH} 2 0.6 ° Experiment iH} 2 0.6 
~3 04 Wind-US SST-V 2 od 

n =, oF ind- ef n 

Zo fates FUN3D SA is 
> 02 — — FUNS3D SST-V => 02 

0.0 0.0 

-100 -50 O 50 100 150 200 250 300 —100 —50 50 100 150 200 250 300 
u velocity, m/s u velocity, m/s 

1.2 xX/H = 6.4 1.2 x/H =7.5 
ae 3x 1.0 
8B 8608 
ue 1 £ 0.6 

=) =) 
iB 22° 
> > 0.2 
0 


“50 


0 50 


100 
u velocity, m/s 


150 200 250 


300 


100 150 200 250 300 
u velocity, m/s 


0. 
—50 0 50 


Figure 4.—Transonic diffuser. Velocity profiles computed with Wind-US and FUN3D. 
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The FUN3D results computed with the Wilcox RSM and the SSG/LRR RSM are plotted with the 
SA and SST-V results in Figure 5 to Figure 7. In the pressure plots of Figure 5, the RSM results are 
slightly better at predicting the pressure rise just downstream of the shock, however, the values further 
downstream are slightly higher than the experimental values, whereas the SA and SST-V results agree 
with the experimental values for x/Htprogt Values between 6.5 to 9. The velocity profiles are shown 
in Figure 6. The RSMs do not predict separation, but they do a better job of predicting the boundary 
layer shape at x/Henroat = 6.4 and 7.5. The RMS u-velocity is shown in Figure 7. It was computed 
by taking the square root of the u'w’ from the FUN3D output. For the SA model, uw’ is computed 
from the Boussinesq approximation (Eq. (1)) using Bradshaw’s approximation for k, where 
k = Ur, 25;;5;;/(0.31p). The plots show that the RSMs are better than the SA and SST-V models at 
predicting the values at the peaks of the profiles. For x/Hip7oq¢ = 2-9 and 4.6, the values at the center of 
the diffuser are under-predicted for all of the models. 


1.0 | 


g | Bottom wall g Top wall 

2p 0.9) 2m 0.9) 

o 2 oo ° 

50 0.8 50 0.84 

_—o =o 

wr 07 a- 07] 

= ll = ll 

Be 0.6 | 3 & 0.64 

= Pu 0.5 a Pa 0.5 

ES 04 — Wilcox RSM ES o4 

fe) - ++ SSG/LRR RSM {e) 

2 0.3 z 0.3 ; : : 
- -3 -1 1 3 5 i 9 - -3 -1 1 3 5 7 9 

X/H throat: throat = 0.1445 ft X/Hthroat throat = 0.1445 ft 


Figure 5.—Transonic diffuser. Normalized wall pressure on the bottom and top walls computed with FUN3D. 
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Figure 6.—Transonic diffuser. Velocity profiles computed with FUN3D. 
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Figure 7.—Transonic diffuser. RMS u-velocity computed with FUN3D. 


Overall, it is unclear if the RSMs offer benefits over the traditional turbulence models based solely on 
this case. For the pressure, they are better at predicting the pressure rise just downstream of the shock, but 
then over-predict it further downstream. For the velocity profiles, the RSMs are worse than the traditional 
models at predicting the separation on the top wall, but give better agreement with data further 
downstream. For the u,;.,,; profiles, the RSMs are significantly better at predicting the profiles. 


4.2 Axisymmetric Supersonic Compression Corner 


4.2.1 Experimental Setup 


This study examines supersonic flow over an axisymmetric compression corner (Refs. 11, 13, 15, 27, 
and 28). The model is a 5.08 cm diameter cylinder with a 30-degree flare, which generates a shock wave 
(Figure 8). The cylinder has an upstream cusped nose designed to minimize the strength of the shocks, 
and data indicates that reflected shocks from the tunnel walls have no effect in the measurement region of 
interest. The flare is located 1.0 m downstream of the cusp-tip, allowing a turbulent boundary layer to 
develop upstream of the shock-wave boundary-layer interaction. The flare surface begins at x = 0.0 cm 
and ends at x = 5.196 cm. The test section has a Mach number of 2.85, a Reynolds number of 16x10°/m, a 
stagnation pressure of 1.7 atm and a stagnation temperature of 270 K. The data consists of laser Doppler 
velocimeter mean velocities and Reynolds stresses, surface static pressures, Schlieren photography, oil 
flow visualizations, and holographic interferometry data. 

This flow has a complex shock-wave/boundary-layer interaction region, as shown by the FUN3D 
results in Figure 9, depicting Mach contours in the interaction region computed with the SA model. A 
primary shock wave is generated when the flow encounters the flare surface. The shock penetrates the 
subsonic portion of the boundary layer causing a strong adverse pressure-gradient, which leads to 
separation in the corner between the cylinder and the flare surfaces. The separation extends upstream of 
the flare, inducing a secondary shock. On the flare surface, the flow at the reattachment point is turned 
parallel to the flare, generating a shock that coalesces with the primary shock. The flow is turned back 
axially downstream at the cone-afterbody, causing the flow to expand and weaken the shock wave. 
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Figure 9.—Axysimmetric compression corner. Mach contours 
computed using the SA turbulence model. 


4.2.2 CFD Methodology 


For this case, a single-cell, axisymmetric wedge grid was used for computations with the SA and 
SST-V models. Since wedge boundary conditions are not yet available in FUN3D for the RSMs, a 
90-degree version of the grid was generated for use with the RSMs. The axisymmetric grid was originally 
generated by Debonis (Ref. 12) using Pointwise (Ref. 29) for structured Wind-US computations. It begins 
on the cylinder body, downstream of the cusp. The length of the cylinder, 75 cm, was chosen such that the 
CFD matched the experimental boundary-layer profile at the first measurement station. The grid was 
thoroughly tested for grid convergence in Reference 12, resulting in a dense grid in order to capture the 
shock and separation features. The grid was orthogonal to the wall, with wall spacing of 5.0x10~° cm and 
a maximum y’ of 0.2. Grid lines are parallel to the shock location with clustering near the expected shock 
location. There are a total of 1,265 axial points and 729 radial points for a total of 922,185 points. 
Pointwise was used to save the axisymmetric grid in unstructured, hexahedral format compatible with 
FUN3D. For the computations using the RSMs, the axisymmetric grid was extruded 90° about the x-axis 
at 5.6° increments for a total of 17 equally-spaced points in the circumferential direction and a total of 
15,478,857 points. The grid is shown in Figure 10. 
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Figure 10.—Axisymmetric compression corner. Grid. (a) Entire domain. (b) Closeup of flare, every 8th point shown. 


The FUN3D boundary conditions were set equivalent to those used by DeBonis for Wind-US. The 
inflow boundary was fixed using the supersonic freestream conditions. The outflow was extrapolated. The 
upper boundary was a farfield condition which used Riemann invariants. No-slip adiabatic viscous 
conditions were set at the solid walls. The Roe inviscid flux method (Ref. 30) with the min-mod inviscid 
flux limiter (Ref. 31) was initially used for all of the cases. FUN3D was first run with the SA and SST-V 
models, and converged results, where solutions were considered to be converged with the most sensitive 
parameters of the solution, the skin friction and the shear stress profiles, unchanging, were compared with 
the Wind-US results of DeBonis. The results from both models were in good agreement between the two 
codes, confirming that the FUN3D inputs and boundary conditions were set to appropriate values. Next 
solutions were obtained with the SSG/LRR and Wilcox RSMs using similar input parameters. These 
solutions required smaller CFL numbers and more iterations to reach convergence. The solution residuals 
were low but more oscillatory than desired. In an attempt to reduce these oscillations, the solution was 
then run using the Van Albada flux limiter with a heuristic pressure limiter (hvanalbada) (Ref. 33), with 
the limiter frozen after 20000 iterations. The resulting solutions had slightly less oscillation of the 
residuals, and the results were in better agreement with the experiment. The SST-V model was also run 
with the same Van Albada flux limiter settings, and the residuals and the solution were comparable to the 
solution run with the minmod limiter. In the results which follow, the RSM results in best agreement with 
the data are shown. It must be emphasized though, that without full iterative convergence, these results 
are incomplete, due to the equations of motion not being entirely satisfied. 


4.2.3 Results 


For this case, we examined the static pressure and skin friction coefficient on the model, axial and 
radial velocity profiles, and shear stress profiles. The surface static pressure is shown in Figure 11(a). The 
first pressure rise upstream of the compression corner is caused by the shock-wave induced corner 
separation. The larger pressure rise beginning at the compression corner (x = 0 cm) is caused by the shock 
wave induced by the flare. The pressure increases until it reaches the end of the flare, then drops as the 
flow expands due to the flow turning back axially. The SA model shows the best agreement with the 
experimental pressure values. It predicts the start of the first pressure rise accurately, and predicts the 
correct pressure rise at the corner and does well on the flare. The small blip in pressure at the corner is a 
result of a stagnation region caused by a counter-rotating vortex in the compression corner. The SST-V 
model predicts the first pressure rise at approximately x = —4.4 cm, which is upstream of the experimental 
value, and under-predicts the maximum pressure rise on the flare slightly. The SSG/LRR RSM predicts 
the first pressure rise slightly farther downstream than the SA model, at x = —2.6 cm, and predicts an 
overly steep pressure rise on the flare. The Wilcox RSM predicts the onset of separation upstream of the 
experiment at x = —3.6 cm, and the pressure rise on the flare is in good agreement with the experiment. 
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The skin friction coefficient, Cy, is plotted in Figure 11(b) and indicates the attached and separated 
regions in the flow. Upstream of the separation, the reported range of experimental C; values is 0.00155 
to 0.0017, (Ref. 28), as indicated by the black lines in Figure 11(b), and the separation begins between 
—3.25 and —2.75 cm (Ref. 12). All turbulence models predict the correct skin friction coefficient values 
upstream of the separation. The separation locations correspond closely with the first pressure rise of 
Figure 11(a). The SA model solution shows the separation beginning at x = —2.74 cm, which agrees with 
the experiment. The SST-V model and Wilcox RSM predict the onset of separation upstream of the 
experimental value, at -4.40 and —3.58 cm, respectively. The SSG/LRR RSM predicts beginning of the 
separation downstream of the experiment at —2.34 cm. All of the models show a rise in Cy, in the corner 
due to the counter-rotating vortex. The Wilcox and SSG/LRR RSMs reattach at x = 1.26 and 0.52 cm, 
respectively, and the SA and SST-V models reattach further downstream at 2.33 and 2.44 cm. The spike 
in Cy, at the aft end of the flare, is due to the high shear stress as the flow accelerates and expands around 
the corer. The separation and reattachment values are summarized in Table 3. None of the models do 
well at predicting the separation length, however both RSMs predict separation lengths closer to the 
experimental value than the SA and SST-V models, with the SSG/LRR RSM predicting the best value, 
22.7 percent smaller than the experiment. 
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Figure 11.—Axisymmetric compression corner. Surface profiles. (a) Static pressure. (b) Skin friction 
coefficient. 


TABLE 3.—AXISYMMETRIC COMPRESSION CORNER. SEPARATION AND REATTACHMENT LOCATIONS 
Separation location, Reattachment location, Separation length, A(Separation length)exp, 
cm cm cm percent 
Experimental oil flow —2.730 0.970 3.700 J wenenn= 
(approximate values) 
SA —2.74 2.33 5.07 37.0 
SST-V —4.40 2.44 6.84 84.9 
Wilcox RSM -3.58 1.26 4.84 30.8 
SSG/LRR RSM —2.34 0.52 2.86 —22.7 
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The axial velocity profiles are shown in Figure 12. At the first measurement station, at x = 4.5 cm, 
the flow is attached and the velocity profile thickness of the computations matches the experiment. All 
turbulence models agreed well with each other and with the data, though they were slightly less full at the 
wall. At x =—3.0 cm, the solutions computed with the SST-V and Wilcox models are separated, and the 
SA model and SSG/LRR RSM solutions and the experimental boundary layers remain attached. All of the 
models show separation at x = —2.0 cm, where the experiment is also separated. The SA, SST-V, and 
Wilcox RSM give a stronger and larger separated region than the experiment, whereas the SSG/LRR 
RSM gives a smaller separated region. Once the flow has separated, there is very little agreement between 
models and none of the models predict the velocity profiles well, with the SST-V model showing the 
strongest and largest separation. However, once the flow attaches and turns the corner to flow in the axial 
direction, all of the profiles approach the experiment, as shown at x = 5.896 cm, with the SSG/LRR result 
still deviating the most. 

The turbulent shear stress profiles are shown in Figure 13. The SA and SST-V results are similar to 
those computed with Wind-US and reported by DeBonis (Ref. 12). At the most upstream station, x = 
—4.5 cm, all of the models under-predict the shear stress near the wall, but agree with the experiment in 
the outer portion of the boundary layer. Overall, the CFD shear stress profiles are in poor agreement with 
the experiment on the flare. Once the flow has turned back axially, at x = 5.896 cm, the profiles produced 
by the RSMs show improvement over the standard turbulence models. 

Overall, for the supersonic axisymmetric compression corner, the RSMs were more difficult to 
converge and more sensitive to code inputs than the standard turbulence models. The SA model does the 
best job of computing the onset of separation. The RSMs did better at predicting the separation length 
than the standard models. None of the models did well at predicting the velocity and shear stress profiles 
in the separated region. For this case, use of the RSMs was found to provide limited benefit. 


4.3 Mixing Layer 
4.3.1 Experimental Configuration 


In this case, two supersonic streams, initially separated by a splitter plate, come into contact and form 
a shear layer. Results are compared with test case 2, in the experiment of Goebel and Dutton (Ref. 14). At 
the entrance of the mixing layer, the flow conditions are M,=1.91, T;,=578 K, and U;=700 m/s for the 
high-speed stream and M,=1.36, T;,=295 K, and Uz=399 m/s for the low-speed stream. Both streams 
have a pressure of 49 kPa. The relative Mach number, M, = AU/(@) = 0.91, where AU = (U, — Uj), 
and a is the average of the freestream speed of sound. 


4.3.2 CED Methodology 


The Wind-US computations shown below are nearly identical to those reported by Yoder (Ref. 20), 
with only some minor differences in the grid. The two-dimensional structured grid consisted of three 
zones. Zones | and 2 were upstream of the mixing layer, each having 101 points in the axial direction and 
89 points in the vertical direction. The lengths of the incoming streams were chosen to match the 
experimental values of the momentum thickness. Lengths of 240 and 160 mm were required for the 
high-speed and low-speed streams, respectively. The third zone, beginning at the end of the splitter plate 
was 401x241 and 500 mm in the axial direction by 48 mm in the vertical direction. The plate was 0.5 mm 
thick and had 64 points along the vertical edge. It is shown in Figure 14. The grid was clustered vertically 
on the splitter plate to capture the boundary layers with wall spacing of 0.0025 mm resulting in y* values 
of 0.8 to 1.0. It was also clustered axially near the splitter plate trailing edge such that x” was 
approximately 16. Results computed using a version of this grid with every other point removed were 
nearly equivalent to results computed using the full grid. Viscous adiabatic wall boundary conditions 
were used on the splitter plate. In the experiment, the top and bottom test section walls were angled to 
account for boundary layer growth; in the computations, inviscid boundary conditions were used. At the 
inflow, the flow was fixed, and at the downstream boundary, the flow was extrapolated. 
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Figure 12.—Axisymmetric compression corner. x-velocity profiles. 
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Figure 13.—Axisymmetric compression corner. Turbulent shear stress profiles. 
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Figure 14.—Mixing layer. Computational mesh. (Every fourth grid point shown.) 


The grid used for FUN3D was created from the Wind-US grid, with a few modifications. It was 
extruded 1 mm in the z-direction, to make it three-dimensional. The tight axial packing at the inflow 
boundary was not compatible with FUN3D, so a short region (10 mm) with inviscid walls was added 
upstream of the viscous walls and the axial clustering was relaxed. The total number of axial grid points 
remained the same. 

For the Wind-US computations, the SST-V model and an Explicit Algebraic Stress Model (EASM) 
(Ref. 20) were used. The EASM is a k-epsilon formulation that uses a nonlinear equation to compute the 
shear stress tensor, and is described in Reference 20. Results computed using this model are included for 
comparison because it gives a more complete calculation of the shear stress tensor compared to models 
based on the Boussinesq approximation while remaining less computationally expensive than the RSMs. 
For the FUN3D computations, the SST-V model, the SSG/LRR, and Wilcox RSMs were used. 


4.3.3 Results 


Mean velocity, turbulent intensity, and turbulent shear stress profiles at several axial stations are 
given in Figure 15 to Figure 18, and the shear layer thickness is given in Figure 19. The Wind-US SST-V 
results are also included in these plots and are in excellent agreement with the FUN3D SST-V results. At 
the entrance to the mixing section, Figure 15 to Figure 18 show that all of the turbulence models give 
nearly the same profile and agree well with data, however downstream, the results deviate. The velocity 
profiles are shown in Figure 15 and indicate that, beginning at x = 100 mm, the SST-V model has a linear 
slope in the mixing region, with a sharp change to the freestream values at both edges. The EASM model 
has a nonlinear transition at the edges of the mixing region, and the Wilcox and SSG/LRR RSMs behave 
nearly the same and their profiles fall somewhat in between the SST-V and EASM profiles for the 
curvature at the mixing-layer edges. The streamwise turbulence intensity is shown in Figure 16. At all 
axial stations, the SST-V model under-predicts the maximum values, while the EASM does well at 
computing the maximum values in the center of the mixing section. The Wilcox model predicts similar 
maximum values, while the SSG/LRR RSM over-predicts the maximum values. The profiles are fairly 
similar in shape, with the EASM model having the most gradual transition at the edges of the mixing 
region. The transverse turbulence intensity is given in Figure 17. Again, at the entrance to the mixing 
section, the models all give essentially the same result. Moving downstream, the SST-V model 
significantly over-predicts the maximum values. The EASM and the RSMs do a good job predicting the 
peak values downstream, with the EASM having a more gradual slope at the edges of the mixing region 
and giving a slightly more narrow mixing region. The turbulent shear stress is shown in Figure 18. The 
EASM does the best job of predicting the maximum values and the shape. The two RSMs give very 
similar results with slightly higher maximum values; the SST-V results are very similar, with a slightly 
wider mixing region. 


NASA/TM—2017-219468 18 


20 


10 


Y- Yo, mm 
oO 


20 


10 


Y- Yo, mm 
(oe) 


x=10 mm} 25 50 100 150 


= : | 
i f ; p oe 
aa San are" eke 0 a wf é ‘ 9} 
oh By %, ‘a4 4 g 
: 4 : : : 
: d 5 
Oo o Goebel 
FUN3D SST-V 
a= s as FUN3D Wilcox RSM 
aii ‘= FUN3D SSG/LRR RSM 
commana Wind-US SST-V 
=== Wind-US EASM 
500 1000 1500 2000 
U, m/s (subsequent profiles shifted by 300 m/s) 
o o 
250 300 350 4 400 ¢ 450 © 
x = 200 mm ; } G q 
‘ : q oC a Fi 
3 1) G J! oF H 
a) A 4 fi é by a 
J ) g Y Pa j 4/ i, 
fl 4 i) ds 4 E j / 
# i i fp of £ 
of o 9 ) i & 
i i , 4 ; 4 : 7 A 
7 i ip ae : 
: i ; i" ai : 
< CJ r) B 
4 i U) iD 0 , 
‘ ) : zt 
GO} o = ° 5 
fa) 
500 1000 1500 2000 


U, m/s (subsequent profiles shifted by 300 m/s) 


Figure 15.—Mixing layer. Mean velocity profiles. 
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Figure 16.—Mixing layer. Streamwise turbulence intensity. 
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Figure 17.—Mixing layer. Transverse turbulence intensity. 
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Figure 18.—Mixing layer. Turbulent shear stress profiles. 
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Figure 19.—Mixing layer. Shear layer thickness. 


The shear layer thickness is given in Figure 19 and is defined as the distance b, between transverse 
locations where & = ti, — 0.1At and & = &, + 0.1A%. The RSMs do the best job at getting the shear 
layer thickness correct. The SST-V model results have the highest values and the EASM gets the 
downstream values too low. This corresponds to the streamwise velocity and turbulence quantity profiles 
in Figure 15 to Figure 18. The width of the shear layer predicted by the SST-V is largest, with the RSMs 
slightly smaller. The RSMs have slightly more curvature at the edges of the shear layer than the SST-V 
model and less than the EASM. Since the definition of b uses the 10%AU criteria, the curvature at the 
edges of the shear layer effects the value. 

While all of the turbulence models examined appear to predict the streamwise velocity profiles well, 
it appears that the EASM and the RSMs have benefits over the two-equation SST-V model. Overall, they 
are closer to the experimental values in predicting the turbulence intensity, turbulent shear stress and 
shear layer thickness. This is most likely due to the fact that these models account for anisotropy effects 
in their calculation of the shear stress tensor. 


4.4  Axisymmetric Jet 


4.4.1. Experimental Configuration 

In this experiment, subsonic, unheated air is passed through an axisymmetric convergent nozzle to 
produce a Mach number of Mjez = Ujet/Ajer = 0.51, exiting into quiescent air. (The acoustic Mach 
number is, Mg = Uje¢/Qyer = 0.5.) This experiment used an acoustic research nozzle referred to as 
ARN2, having an exit diameter of 2 in., at conditions referred to as set point 3 (Ref. 35). The freestream 
static temperature and pressure were 530 °R and 14.3 psi, respectively. Particle image velocimetry (PIV) 
data was taken of axial and turbulent velocities. 


4.4.2 CFD Methodology 


This case was computed using the Wind-US code with the SA and SST-V turbulence models as 
described on the TMR (Ref. 19). Since quiescent air is often difficult to compute with CFD, the flow 
external to the nozzle was set to M,¢¢ = 0.01, where the external flow conditions are referred to as the 
reference conditions. At the jet inflow, p;/Prer = 1.197 and T;/T;yer = 1.0. The nozzle walls had 
adiabatic viscous boundary conditions and the jet outflow was extrapolated. A series of grids was 
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available on the TMR, each coarse grid was created by removing every other grid point from the next 
finest grid. The grids contained three zones each as shown in Figure 20 and Table 4. To convert these 
grids to FUN3D-compatible format, Pointwise was used, resulting in a single zone grid. The grids were 
one-cell deep spanning approximately 1° in the theta direction. For computations using the RSMs, a 90° 
sector was required. To generate this grid, the axisymmetric grid was extruded 90° in the theta direction at 
5-degree increments for a total of 19 grid points. A grid resolution study was done using FUN3D with the 
SA and SST-V turbulence models. The solution did not show appreciable changes when using Grid 129 
and Grid 257, so Grid 129 was used for the remainder of the computations. 

Regardless of the turbulence model chosen, none of the solutions converged to a definite steady state 
result. When run in a time-accurate mode, solutions using the SA, SST-V and SSG/LRR model 
converged to a quasi-steady result when the solutions were time-averaged, however the solution 
computed using the Wilcox RSM was still noticeably changing after over one million time steps, and so 
the results are not included in the plots. The solutions presented on the TMR were computed with the 
CFL3D (Ref. 36) and the Wind-US codes, and also required a time-accurate solution method to achieve a 
quasi-steady state result. Wind-US results using the SA and SST-V turbulence models on Grid 257, were 
also given on the TMR. FUN3D results computed on Grid 257 computed using the SA and SST-V 
models were in good agreement with these Wind-US results. 
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Figure 20.—Axisymmetric jet. Grid shown at the second from the finest level (Grid 129). 


TABLE 4.—AXISYMMETRIC JET 


[Grid study. ] 
Name Upstream external region | Internal nozzle region | Jet mixing region 
Grid 65 25x25 16x25 65x57 
Grid 129 49x49 31x49 129x113 
Grid 257 97x97 61x97 257x225 
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4.4.3 Results 


The results are given in Figure 21 to Figure 26 and include centerline velocity and turbulent kinetic 
energy profiles, and radial velocity, turbulent shear stress and turbulent kinetic energy profiles at five 
locations downstream of the jet exit. As mentioned above, the solution computed using the Wilcox RSM 
did not achieve a quasi-steady result, and the results are omitted from the plots. The mean x-velocity is 
given in Figure 21 and indicates that the SSG/LRR RSM does the best job of predicting the overall 
mixing. Although the onset of mixing takes effect at about 6.5 diameters downstream of the jet exit, 

2 nozzle-exit diameters downstream of experimental onset location of x/Dje; = 4.5, the SSG/LRR RSM 
result otherwise gives the best overall mixing throughout the mixing region. The SA model result shows 
the onset of mixing closest to the experiment at x/Dje¢ = 5.8, however it over-mixes, resulting in a lower 
downstream centerline velocity. The SST-V model result shows the onset beginning the furthest 
downstream of the experiment at x/Dje, = 8.1, then over-mixes, resulting in downstream velocities 
similar to the SA result. The turbulent kinetic energy, where k for the SA result was computed using the 
Boussinesq approximation and Bradshaw’s approximation as described in Section 4.1.3, is plotted in 
Figure 22. The general profile shape is similar to that of the experiment. However, the values begin to 
increase further downstream, and the peak values are quite different. The SA model peak turbulent kinetic 
energy is 60 percent lower than the experiment, the SST-V peak value is 22 percent higher and the 
SSG/LRR peak is 42 percent higher. All of the CFD results show that the turbulent kinetic energy begins 
to increase between 5 and 8 jet diameters downstream of the jet exit, whereas the experimental turbulent 
kinetic energy begins to rise at the jet exit. For the SA model, the slope as the turbulent kinetic energy 
increases is closest to the experiment; the SST-V and SSG/LRR RSM have much steeper slopes. 
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Figure 21.—Axisymmetric jet. Axial velocity profiles on 
the centerline computed using FUN3D. 
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Figure 22.—Axisymmetric jet. Turbulent kinetic energy profiles 
along the centerline computed using FUN3D. 
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Figure 23.—Axisymmetric jet. Radial x-velocity profiles computed using 
FUNS3D. (Subsequent profiles shifted by u/Uje¢ = 1.0.) 


The x-velocity profiles are shown in Figure 23 at x/Dje¢ = 2.,5.,10.,15.,and 20. Near the jet exit, 
all of the profiles agree well with the experiment and predict the correct profile shapes. Beginning 
at x/Dje¢ = 10., the values near the centerline are less than the experiment for the SA and SST-V 
computations, as was shown in Figure 21, but agree well away from the centerline. The y-velocity, 
turbulent shear stress and turbulent kinetic energy are shown in Figure 24 to Figure 26, and show how the 
computations predict the shear layer mixing. At x/Dje, = 2. and 5., the SSG/LRR RSM does slightly 
better at predicting the profile shape. However, further downstream, the y-velocity and the turbulent shear 
stress results are similar for all three models. At the three most downstream stations, the turbulent kinetic 
energy is different at the centerline, as was also indicated by the centerline profile in Figure 22. Overall 
for this subsonic axisymmetric jet study, the SSG/LRR model shows some benefits over the SA and 
SST-V models at predicting the mixing. 
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Figure 24.—Axisymmetric jet. Radial y-velocity profiles computed using FUN3D. 
(Subsequent profiles shifted by v/Uj., = 0.04.) 
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Figure 25.—Axisymmetric jet. Radial turbulent shear stress profiles computed 
using FUN3D. (Subsequent profiles shifted by u'v'/Uje:” = 0.01.) 
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Figure 26.—Axisymmetric jet. Radial turbulent kinetic energy profiles computed 
using FUN3D. (Subsequent profiles shifted by k/Ujet" = 0.03.) 


5.0 Conclusion 


Two full second-moment Reynolds stress turbulence models available in the FUN3D code, the 
Wilcox and the SSG/LRR, were evaluated for four test cases: a transonic diffuser, a supersonic 
axisymmetric compression corner, a supersonic compressible planar mixing layer, and a subsonic 
axisymmetric jet. These model results were compared with solutions computed using the SA and SST-V 
turbulence models, and for the planar mixing layer case, an EASM. The transonic diffuser was the first 
case evaluated, and gave inconclusive results as to the benefits of the RSMs. The supersonic 
axisymmetric compression corner, which has a complex shock system and a large region of separated 
flow, was evaluated next. The SA model did the best at computing the pressure rise and the onset of 
separation. The RSMs were better at predicting the separation length. All of the models had difficulty 
computing the boundary layer profiles and turbulence quantities in the separated region. It was more 
difficult to obtain converged solutions using the RSMs and solutions required more iterations, therefore 
limited benefit was gained by using the RSMs for this case. For the supersonic planar mixing layer, both 
RSMs and the EASM had benefits over the SST-V model at predicting the turbulence intensity, turbulent 
shear stress and shear layer thickness. These benefits were most likely because these models account for 
anisotropy effects in their calculation of the shear stress tensor. For the subsonic axisymmetric jet, the 
SSG/LRR model predicted the mixing of the core velocity the best, whereas the solution with the Wilcox 
RSM did not reach a stable answer for this case. 

The four cases examined are flows which are challenging for current turbulence models because they 
contain mixing, shock waves and/or separation. Overall, the RSMs showed benefit over the SA and 
SST-V models for the planar mixing layer and the axisymmetric jet flow, and may be useful for future 
nozzle calculations. The EASM, which is less cumbersome to code and requires less computational time, 
did well with the planar mixing layer, and also may be a viable choice for mixing flows. While the cases 
examined are challenging flows, they are still relatively simple in geometry and flow features. More 
complex flow cases may reveal more benefits of the RSMs and are recommended for future study. 
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